rm(list = ls())
library(foreign)
library(splm)
library(spdep)
library(spam)
library(ape)
library(dplyr)
library(cshapes)
load("data/data_imputed.RData")


data  <- data_imputed
data <- data %>% 
       dplyr::rename(country = country_name) 

data <- data %>% 
       dplyr::arrange(cowcode, year) %>% 
       dplyr::group_by(cowcode) %>% 
       dplyr::mutate(laglogidealdistance =  dplyr::lag(logidealdistance, n = 1L),
                     lagagree = dplyr::lag(agree, n=1L))

data <- data %>% 
       dplyr::filter(!is.na(laglogidealdistance))

load("data/state_year.RData")


years <- unique(data$year)

moran2 <- list()
for (i in seq_along(years)) {
       
       cowcodes <- data %>% dplyr::filter(year == years[i]) %>% 
              dplyr::pull(cowcode)
       df_slice <- data %>% dplyr::filter(year == years[i])
       slice <- state_year[state_year$year == years[i],]
       slice <- slice[slice$cowcode %in% cowcodes,]
       
       slice <- as_Spatial(slice)
       coords <- coordinates(slice)
       kn5 <- knn2nb(knearneigh(coords, k = 5))
      listwk5 <- nb2listw(kn5, style="W", zero.policy=TRUE)
      W_k5 <-  as.matrix(as(as_dgRMatrix_listw(listwk5), "CsparseMatrix"))
      moran <- Moran.I(df_slice$laglogidealdistance, weight = W_k5)
      moran2[[i]] <- moran
       
}


# Generate matrix of results
mat <- matrix(NA, nrow = 4, ncol=length(years))
colnames(mat) <- years
for (i in seq_along(moran2)){

       mat[,i] <- unlist(moran2[[i]])
       rownames(mat) <- c("Observed", "Expected", "sd","p.value")
}

mat <- t(mat)
mat_df <- as.data.frame(mat)
mat_df$year <- rownames(mat)

library(ggplot2)

library(tid)
df_moran <- gather(mat_df, key = "Type", value = "Moran", Observed:Expected)

ggplot(df_moran, aes(x=year, y= Moran, group = Type, linetype = Type)) +
       geom_point() + 
       geom_line() +
       scale_x_discrete(breaks = seq(1979, 2014, 3)) +
       labs(x = "Year", y = "Moran's I") +  
       theme_bw() + 
       theme(legend.title =  element_blank(),
             legend.position = "bottom",
             text = element_text(size=16))
       
ggsave("moranI22.pdf", height =5 , width = 10)


####Local # Local Moran's I
listw <- nb2listw(cshp_10_nbBin_10$neighbours, style="W") 
s <- list()
for (i in 1:length(years)){
s[[i]] <- localmoran(data_subset[[i]]$s2dfidealpoint,listw,p.adjust.method="bonferroni") %>% 
            as.data.frame() %>%
              dplyr::rename(p.value = `Pr(z > 0)`) %>% 
            dplyr::filter(p.value < 0.1) %>% dplyr::select(Z.Ii) %>% dplyr::mutate(year = years[i]) %>% 
            data.frame()
}
local <- list()
for (i in seq_along(years)) {
       
       cowcodes <- data %>% dplyr::filter(year == years[i]) %>% 
              dplyr::pull(cowcode)
       df_slice <- data %>% dplyr::filter(year == years[i])
       slice <- state_year[state_year$year == years[i],]
       slice <- slice[slice$cowcode %in% cowcodes,]
       
       slice <- as_Spatial(slice)
       coords <- coordinates(slice)
       kn5 <- knn2nb(knearneigh(coords, k = 5))
       listwk5 <- nb2listw(kn5, style="W", zero.policy=TRUE)
       loc_moran <- spdep::localmoran(x = df_slice$laglogidealdistance, listw = listwk5) %>% 
           as.data.frame() %>%
           dplyr::rename(p.value = `Pr(z != E(Ii))`) %>% 
           dplyr::filter(p.value < 0.1) %>% 
           dplyr::select(Z.Ii) %>% dplyr::mutate(year = years[i]) 
       local[[i]] <- loc_moran
       
}

##################################################
# Histograms of significant local Morans I
df_localmoran <- do.call(rbind, local)
ggplot(df_localmoran , aes(factor(year), Z.Ii)) + 
       geom_boxplot(aes(fill = factor(year))) + labs(x = "Year", y = "Z-scores") + 
       scale_x_discrete(breaks = seq(1979, 2014, 3)) +
       theme_bw() + theme(panel.background = element_blank(), panel.grid.major = element_blank(), 
                          panel.grid.minor = element_blank(), axis.ticks.length = unit(0,"cm"),
                          legend.position = "none",
                          axis.text=element_text(size=10)) 
ggsave("localmoran2.pdf", height =5 , width = 10)



